pydaptivefiltering

134class AdaptiveFilter(ABC):
135    """Abstract base class for all adaptive filters.
136
137    Parameters
138    ----------
139    filter_order:
140        Order in the FIR sense (number of taps - 1). For non-FIR structures, it can be used
141        as a generic size indicator for base allocation.
142    w_init:
143        Initial coefficient vector. If None, initialized to zeros.
144
145    Notes
146    -----
147    - Subclasses should set `supports_complex = True` if they support complex-valued data.
148    - Subclasses are expected to call `_record_history()` every iteration (or use helper methods)
149      if they want coefficient trajectories.
150    """
151
152    supports_complex: bool = False
153
154    def __init__(self, filter_order: int, w_init: Optional[ArrayLike] = None) -> None:
155        self.filter_order: int = int(filter_order)
156        self._dtype = complex if self.supports_complex else float
157
158        self.regressor: np.ndarray = np.zeros(self.filter_order + 1, dtype=self._dtype)
159
160        if w_init is not None:
161            self.w: np.ndarray = np.asarray(w_init, dtype=self._dtype)
162        else:
163            self.w = np.zeros(self.filter_order + 1, dtype=self._dtype)
164
165        self.w_history: List[np.ndarray] = []
166        self._record_history()
167
168    def _record_history(self) -> None:
169        """Store a snapshot of current coefficients."""
170        self.w_history.append(np.asarray(self.w).copy())
171
172    def _final_coeffs(self, coefficients: Any) -> Any:
173        """Return last coefficients from a history container (list or 2D array)."""
174        if coefficients is None:
175            return None
176        if isinstance(coefficients, list) and len(coefficients) > 0:
177            return coefficients[-1]
178        try:
179            a = np.asarray(coefficients)
180            if a.ndim == 2:
181                return a[-1, :]
182        except Exception:
183            pass
184        return coefficients
185
186    def _pack_results(
187        self,
188        outputs: np.ndarray,
189        errors: np.ndarray,
190        runtime_s: float,
191        error_type: str = "a_priori",
192        extra: Optional[Dict[str, Any]] = None,
193    ) -> OptimizationResult:
194        """Centralized output packaging to standardize results."""
195        return OptimizationResult(
196            outputs=np.asarray(outputs),
197            errors=np.asarray(errors),
198            coefficients=np.asarray(self.w_history),
199            algorithm=self.__class__.__name__,
200            runtime_ms=float(runtime_s) * 1000.0,
201            error_type=str(error_type),
202            extra=extra,
203        )
204
205    def filter_signal(self, input_signal: ArrayLike) -> np.ndarray:
206        """Filter an input signal using current coefficients.
207
208        Default implementation assumes an FIR structure with taps `self.w` and
209        regressor convention:
210            x_k = [x[k], x[k-1], ..., x[k-m]]
211        and output:
212            y[k] = w^H x_k   (Hermitian for complex)
213        """
214        x = np.asarray(input_signal, dtype=self._dtype)
215        n_samples = x.size
216        y = np.zeros(n_samples, dtype=self._dtype)
217
218        x_padded = np.zeros(n_samples + self.filter_order, dtype=self._dtype)
219        x_padded[self.filter_order:] = x
220
221        for k in range(n_samples):
222            x_k = x_padded[k : k + self.filter_order + 1][::-1]
223            y[k] = np.dot(self.w.conj(), x_k)
224
225        return y
226
227    @classmethod
228    def default_test_init_kwargs(cls, order: int) -> dict:
229        """Override in subclasses to provide init kwargs for standardized tests."""
230        return {}
231
232    @abstractmethod
233    def optimize(
234        self,
235        input_signal: ArrayLike,
236        desired_signal: ArrayLike,
237        **kwargs: Any,
238    ) -> Any:
239        """Run the adaptation procedure.
240
241        Subclasses should return either:
242        - OptimizationResult (recommended), or
243        - dict-like with standardized keys, if you are migrating older code.
244        """
245        raise NotImplementedError
246
247    def reset_filter(self, w_new: Optional[ArrayLike] = None) -> None:
248        """Reset coefficients and history."""
249        if w_new is not None:
250            self.w = np.asarray(w_new, dtype=self._dtype)
251        else:
252            self.w = np.zeros(self.filter_order + 1, dtype=self._dtype)
253        self.w_history = []
254        self._record_history()

 28class LMS(AdaptiveFilter):
 29    """
 30    Complex LMS (Least-Mean Squares).
 31
 32    Implements the complex LMS recursion (Algorithm 3.2 - Diniz) for adaptive FIR
 33    filtering.
 34
 35    Notes
 36    -----
 37    - Complex-valued implementation (supports_complex = True).
 38    - Uses the unified base API via `@validate_input`:
 39        * optimize(input_signal=..., desired_signal=...)
 40        * optimize(x=..., d=...)
 41        * optimize(x, d)
 42    - This implementation returns the a priori error: e[k] = d[k] - y[k] where
 43      y[k] = w[k]^H x_k, and then updates:
 44          w[k+1] = w[k] + mu * conj(e[k]) * x_k
 45    """
 46
 47    supports_complex: bool = True
 48
 49    step_size: float
 50
 51    def __init__(
 52        self,
 53        filter_order: int,
 54        step_size: float = 1e-2,
 55        w_init: Optional[ArrayLike] = None,
 56    ) -> None:
 57        """
 58        Parameters
 59        ----------
 60        filter_order:
 61            FIR order M (number of taps is M+1).
 62        step_size:
 63            Step-size (mu).
 64        w_init:
 65            Optional initial coefficients (length M+1). If None, zeros.
 66        """
 67        super().__init__(filter_order=int(filter_order), w_init=w_init)
 68        self.step_size = float(step_size)
 69
 70    @validate_input
 71    def optimize(
 72        self,
 73        input_signal: np.ndarray,
 74        desired_signal: np.ndarray,
 75        verbose: bool = False,
 76    ) -> OptimizationResult:
 77        """
 78        Run LMS adaptation.
 79
 80        Parameters
 81        ----------
 82        input_signal:
 83            Input signal x[k].
 84        desired_signal:
 85            Desired signal d[k].
 86        verbose:
 87            If True, prints runtime.
 88
 89        Returns
 90        -------
 91        OptimizationResult
 92            outputs:
 93                Filter output y[k].
 94            errors:
 95                A priori error e[k] = d[k] - y[k].
 96            coefficients:
 97                History of coefficients stored in the base class.
 98            error_type:
 99                "a_priori".
100        """
101        tic: float = perf_counter()
102
103        x: np.ndarray = np.asarray(input_signal, dtype=complex).ravel()
104        d: np.ndarray = np.asarray(desired_signal, dtype=complex).ravel()
105
106        n_samples: int = int(x.size)
107        m: int = int(self.filter_order)
108
109        outputs: np.ndarray = np.zeros(n_samples, dtype=complex)
110        errors: np.ndarray = np.zeros(n_samples, dtype=complex)
111
112        x_padded: np.ndarray = np.zeros(n_samples + m, dtype=complex)
113        x_padded[m:] = x
114
115        for k in range(n_samples):
116            x_k: np.ndarray = x_padded[k : k + m + 1][::-1]
117
118            y_k: complex = complex(np.vdot(self.w, x_k)) 
119            outputs[k] = y_k
120    
121            e_k: complex = d[k] - y_k
122            errors[k] = e_k
123
124            self.w = self.w + self.step_size * np.conj(e_k) * x_k
125
126            self._record_history()
127
128        runtime_s: float = perf_counter() - tic
129        if verbose:
130            print(f"[LMS] Completed in {runtime_s * 1000:.03f} ms")
131
132        return self._pack_results(
133            outputs=outputs,
134            errors=errors,
135            runtime_s=runtime_s,
136            error_type="a_priori",
137        )

 28class NLMS(AdaptiveFilter):
 29    """
 30    Complex NLMS (Normalized Least-Mean Squares).
 31
 32    Implements the normalized LMS recursion (Algorithm 4.3 - Diniz) for adaptive
 33    FIR filtering. The update uses a step normalized by the regressor energy:
 34
 35        y[k] = w[k]^H x_k
 36        e[k] = d[k] - y[k]
 37        mu_k = mu / (||x_k||^2 + gamma)
 38        w[k+1] = w[k] + mu_k * conj(e[k]) * x_k
 39
 40    Notes
 41    -----
 42    - Complex-valued implementation (supports_complex = True).
 43    - Uses the unified base API via `@validate_input`.
 44    - `gamma` is a regularization constant to avoid division by zero.
 45    """
 46
 47    supports_complex: bool = True
 48
 49    step_size: float
 50    gamma: float
 51
 52    def __init__(
 53        self,
 54        filter_order: int,
 55        step_size: float = 1e-2,
 56        gamma: float = 1e-6,
 57        w_init: Optional[ArrayLike] = None,
 58    ) -> None:
 59        """
 60        Parameters
 61        ----------
 62        filter_order:
 63            FIR order M (number of taps is M+1).
 64        step_size:
 65            Base step-size (mu).
 66        gamma:
 67            Regularization to avoid division by (near) zero.
 68        w_init:
 69            Optional initial coefficients (length M+1). If None, zeros.
 70        """
 71        super().__init__(filter_order=int(filter_order), w_init=w_init)
 72        self.step_size = float(step_size)
 73        self.gamma = float(gamma)
 74
 75    @validate_input
 76    def optimize(
 77        self,
 78        input_signal: np.ndarray,
 79        desired_signal: np.ndarray,
 80        verbose: bool = False,
 81    ) -> OptimizationResult:
 82        """
 83        Run NLMS adaptation.
 84
 85        Parameters
 86        ----------
 87        input_signal:
 88            Input signal x[k].
 89        desired_signal:
 90            Desired signal d[k].
 91        verbose:
 92            If True, prints runtime.
 93
 94        Returns
 95        -------
 96        OptimizationResult
 97            outputs:
 98                Filter output y[k].
 99            errors:
100                A priori error e[k] = d[k] - y[k].
101            coefficients:
102                History of coefficients stored in the base class.
103            error_type:
104                "a_priori".
105        """
106        tic: float = perf_counter()
107
108        x: np.ndarray = np.asarray(input_signal, dtype=complex).ravel()
109        d: np.ndarray = np.asarray(desired_signal, dtype=complex).ravel()
110
111        n_samples: int = int(x.size)
112        m: int = int(self.filter_order)
113
114        outputs: np.ndarray = np.zeros(n_samples, dtype=complex)
115        errors: np.ndarray = np.zeros(n_samples, dtype=complex)
116
117        x_padded: np.ndarray = np.zeros(n_samples + m, dtype=complex)
118        x_padded[m:] = x
119
120        for k in range(n_samples):
121            x_k: np.ndarray = x_padded[k : k + m + 1][::-1]
122
123            y_k: complex = complex(np.vdot(self.w, x_k))
124            outputs[k] = y_k
125
126            e_k: complex = d[k] - y_k
127            errors[k] = e_k
128
129            norm_xk: float = float(np.vdot(x_k, x_k).real)
130            mu_k: float = self.step_size / (norm_xk + self.gamma)
131
132            self.w = self.w + mu_k * np.conj(e_k) * x_k
133
134            self._record_history()
135
136        runtime_s: float = perf_counter() - tic
137        if verbose:
138            print(f"[NLMS] Completed in {runtime_s * 1000:.03f} ms")
139
140        return self._pack_results(
141            outputs=outputs,
142            errors=errors,
143            runtime_s=runtime_s,
144            error_type="a_priori",
145        )

 29class AffineProjection(AdaptiveFilter):
 30    """
 31    Complex Affine Projection Algorithm (APA).
 32
 33    Implements Algorithm 4.6 (Diniz) using an affine-projection update with data reuse.
 34
 35    Notes
 36    -----
 37    - This implementation supports complex-valued data (supports_complex=True).
 38    - The base decorator `@validate_input` allows calling optimize with:
 39        * optimize(input_signal=..., desired_signal=...)
 40        * optimize(x=..., d=...)
 41        * optimize(x, d)
 42    """
 43
 44    supports_complex: bool = True
 45
 46    step_size: float
 47    gamma: float
 48    memory_length: int
 49
 50    def __init__(
 51        self,
 52        filter_order: int,
 53        step_size: float = 1e-2,
 54        gamma: float = 1e-6,
 55        L: int = 2,
 56        w_init: Optional[ArrayLike] = None,
 57    ) -> None:
 58        """
 59        Parameters
 60        ----------
 61        filter_order:
 62            FIR order M (number of taps is M+1).
 63        step_size:
 64            Step-size / relaxation factor (mu).
 65        gamma:
 66            Diagonal loading regularization to ensure invertibility.
 67        L:
 68            Data reuse factor (projection order). Uses L+1 past regressors.
 69        w_init:
 70            Optional initial weights (length M+1). If None, initializes to zeros.
 71        """
 72        super().__init__(filter_order=int(filter_order), w_init=w_init)
 73        self.step_size = float(step_size)
 74        self.gamma = float(gamma)
 75        self.memory_length = int(L)
 76
 77    @validate_input
 78    def optimize(
 79        self,
 80        input_signal: np.ndarray,
 81        desired_signal: np.ndarray,
 82        verbose: bool = False,
 83        return_internal_states: bool = False,
 84    ) -> OptimizationResult:
 85        """
 86        Run APA adaptation.
 87
 88        Parameters
 89        ----------
 90        input_signal:
 91            Input signal x[k].
 92        desired_signal:
 93            Desired signal d[k].
 94        verbose:
 95            If True, prints runtime.
 96        return_internal_states:
 97            If True, returns the last regressor matrix and last correlation matrix in result.extra.
 98
 99        Returns
100        -------
101        OptimizationResult
102            outputs:
103                Filter output y[k] (a priori).
104            errors:
105                Error e[k] = d[k] - y[k] (a priori).
106            coefficients:
107                Coefficient history (self.w_history) as a 2D array.
108            error_type:
109                "a_priori".
110            extra (optional):
111                last_regressor_matrix, last_correlation_matrix.
112        """
113        tic: float = perf_counter()
114
115        x: np.ndarray = np.asarray(input_signal).ravel()
116        d: np.ndarray = np.asarray(desired_signal).ravel()
117
118        dtype = complex
119        x = np.asarray(input_signal, dtype=dtype).ravel()
120        d = np.asarray(desired_signal, dtype=dtype).ravel()
121        
122        n_samples: int = int(x.size)
123        m: int = int(self.filter_order)
124        L: int = int(self.memory_length)
125
126        outputs: np.ndarray = np.zeros(n_samples, dtype=dtype)
127        errors: np.ndarray = np.zeros(n_samples, dtype=dtype)
128
129        x_padded: np.ndarray = np.zeros(n_samples + m, dtype=dtype)
130        x_padded[m:] = x
131
132        X_matrix: np.ndarray = np.zeros((L + 1, m + 1), dtype=dtype)
133        D_vector: np.ndarray = np.zeros(L + 1, dtype=dtype)
134
135        last_corr: Optional[np.ndarray] = None
136
137        eye_L: np.ndarray = np.eye(L + 1, dtype=dtype)
138
139        for k in range(n_samples):
140            X_matrix[1:] = X_matrix[:-1]
141            X_matrix[0] = x_padded[k : k + m + 1][::-1]
142
143            D_vector[1:] = D_vector[:-1]
144            D_vector[0] = d[k]
145
146            Y_vector: np.ndarray = X_matrix @ self.w.conj()
147            E_vector: np.ndarray = D_vector - Y_vector
148
149            outputs[k] = Y_vector[0]
150            errors[k] = E_vector[0]
151
152            corr_matrix: np.ndarray = (X_matrix @ X_matrix.conj().T) + (self.gamma * eye_L)
153            last_corr = corr_matrix
154
155            try:
156                u: np.ndarray = np.linalg.solve(corr_matrix, E_vector)
157            except np.linalg.LinAlgError:
158                u = np.linalg.pinv(corr_matrix) @ E_vector
159
160            self.w = self.w + self.step_size * (X_matrix.conj().T @ u)
161            self._record_history()
162
163        runtime_s: float = perf_counter() - tic
164        if verbose:
165            print(f"[AffineProjection] Completed in {runtime_s * 1000:.02f} ms")
166
167        extra = None
168        if return_internal_states:
169            extra = {
170                "last_regressor_matrix": X_matrix.copy(),
171                "last_correlation_matrix": None if last_corr is None else last_corr.copy(),
172            }
173
174        return self._pack_results(
175            outputs=outputs,
176            errors=errors,
177            runtime_s=runtime_s,
178            error_type="a_priori",
179            extra=extra,
180        )

 28class SignData(AdaptiveFilter):
 29    """
 30    Sign-Data LMS (complex-valued).
 31
 32    This is a low-complexity LMS variant where the input regressor is replaced by
 33    its elementwise sign:
 34
 35        y[k] = w^H x_k
 36        e[k] = d[k] - y[k]
 37        w <- w + 2 * mu * conj(e[k]) * sign(x_k)
 38
 39    Notes
 40    -----
 41    - Complex-valued implementation (supports_complex=True).
 42    - Uses the unified base API via `@validate_input`.
 43    - Returns a priori error by default (e[k] computed before update).
 44    """
 45
 46    supports_complex: bool = True
 47
 48    def __init__(
 49        self,
 50        filter_order: int,
 51        step_size: float = 1e-2,
 52        w_init: Optional[ArrayLike] = None,
 53    ) -> None:
 54        """
 55        Parameters
 56        ----------
 57        filter_order:
 58            FIR order M (number of taps is M+1).
 59        step_size:
 60            Step-size (mu).
 61        w_init:
 62            Optional initial coefficients (length M+1). If None, zeros.
 63        """
 64        super().__init__(filter_order=int(filter_order), w_init=w_init)
 65        self.step_size = float(step_size)
 66
 67    @validate_input
 68    def optimize(
 69        self,
 70        input_signal: np.ndarray,
 71        desired_signal: np.ndarray,
 72        verbose: bool = False,
 73        return_internal_states: bool = False,
 74    ) -> OptimizationResult:
 75        """
 76        Run Sign-Data LMS adaptation.
 77
 78        Parameters
 79        ----------
 80        input_signal:
 81            Input signal x[k].
 82        desired_signal:
 83            Desired signal d[k].
 84        verbose:
 85            If True, prints runtime.
 86        return_internal_states:
 87            If True, returns the last regressor sign vector in result.extra.
 88
 89        Returns
 90        -------
 91        OptimizationResult
 92            outputs:
 93                Filter output y[k].
 94            errors:
 95                A priori error e[k] = d[k] - y[k].
 96            coefficients:
 97                Coefficient history stored in the base class.
 98            error_type:
 99                "a_priori".
100        """
101        t0 = perf_counter()
102
103        x = np.asarray(input_signal, dtype=complex).ravel()
104        d = np.asarray(desired_signal, dtype=complex).ravel()
105
106        n_samples = int(x.size)
107        m = int(self.filter_order)
108
109        outputs = np.zeros(n_samples, dtype=complex)
110        errors = np.zeros(n_samples, dtype=complex)
111
112        x_padded = np.zeros(n_samples + m, dtype=complex)
113        x_padded[m:] = x
114
115        last_sign_xk: Optional[np.ndarray] = None
116
117        for k in range(n_samples):
118            x_k = x_padded[k : k + m + 1][::-1]
119
120            y_k = complex(np.vdot(self.w, x_k))
121            outputs[k] = y_k
122
123            e_k = d[k] - y_k
124            errors[k] = e_k
125
126            sign_xk = np.sign(x_k)
127            last_sign_xk = sign_xk
128
129            self.w = self.w + (2.0 * self.step_size) * np.conj(e_k) * sign_xk
130            self._record_history()
131
132        runtime_s = float(perf_counter() - t0)
133        if verbose:
134            print(f"[SignData] Completed in {runtime_s * 1000:.03f} ms")
135
136        extra: Optional[Dict[str, Any]] = None
137        if return_internal_states:
138            extra = {"last_sign_regressor": None if last_sign_xk is None else last_sign_xk.copy()}
139
140        return self._pack_results(
141            outputs=outputs,
142            errors=errors,
143            runtime_s=runtime_s,
144            error_type="a_priori",
145            extra=extra,
146        )